Fluorescence/plastid dynamics
Scatterplot of pigmentation
Use a subset of datapoints
numcells2use <- 200
# Concatenate data
FCDP_name <- paste(CSV_folder,"/FCDP_all_",expt.day[1],".csv",sep = "")
FCDP <- read.csv(FCDP_name)
FCDP.Meso.00 <- FCDP[FCDP$Class=='Mesodinium' & FCDP$Prey=='00',]
if(dim(FCDP.Meso.00)[1] > numcells2use){
rows2use <- sample(c(1:dim(FCDP.Meso.00)[1]),numcells2use)
FCDP.Meso.00 <- FCDP.Meso.00[rows2use,]
}
FCDP.Meso.HP <- FCDP[FCDP$Class=='Mesodinium' & FCDP$Prey=='HP',]
if(dim(FCDP.Meso.HP)[1] > numcells2use){
rows2use <- sample(c(1:dim(FCDP.Meso.HP)[1]),numcells2use)
FCDP.Meso.HP <- FCDP.Meso.HP[rows2use,]
}
for(i in 2:length(expt.day)){
FCDP_name <- paste(CSV_folder,"/FCDP_all_",expt.day[i],".csv",sep = "")
FCDP <- read.csv(FCDP_name)
FCDP.Meso.00.hold <- FCDP[FCDP$Class=='Mesodinium' & FCDP$Prey=='00',]
if(dim(FCDP.Meso.00.hold)[1] > numcells2use){
rows2use <- sample(c(1:dim(FCDP.Meso.00.hold)[1]),numcells2use)
FCDP.Meso.00.hold <- FCDP.Meso.00.hold[rows2use,]
}
FCDP.Meso.HP.hold <- FCDP[FCDP$Class=='Mesodinium' & FCDP$Prey=='HP',]
if(dim(FCDP.Meso.HP.hold)[1] > numcells2use){
rows2use <- sample(c(1:dim(FCDP.Meso.HP.hold)[1]),numcells2use)
FCDP.Meso.HP.hold <- FCDP.Meso.HP.hold[rows2use,]
}
FCDP.Meso.00 <- rbind(FCDP.Meso.00,FCDP.Meso.00.hold)
FCDP.Meso.HP <- rbind(FCDP.Meso.HP,FCDP.Meso.HP.hold)
}
Extract the fluorescent particles – HP-fed
# Step 1: Partition dataset to just fluorescent particles
dataset.fluo <- FCDP.Meso.HP[FCDP.Meso.HP$Ch1.Peak>0,]
dataset.fluo <- dataset.fluo[dataset.fluo$Ch2.Peak>0,]
# Step 2: Use flowPeaks package to fit and visualize the number of clusters of points
invisible(capture.output(fp <- suppressMessages(flowPeaks(dataset.fluo[,c("Ch1.Peak",'Ch2.Peak')],tol=.1,h0=.4))))
# The "capture.output" function captures the printed output of the flowPeaks function
# The "invisible" function suppresses that output so it doesn't automatically print.
# The "suppressMessages" function suppresses the messages output of flowPeaks
# Step 3: Extract cluster assignments
dataset.fluo$peak <- fp$peaks.cluster
# Step 4: Assign clusters to colors (or unknown)
fp.data <- as.data.frame(summary(fp))
fp.data$plastid.color <- 'red' # Presume plastids are red in color
for(i in 1:dim(fp.data)[1]){
if(fp.data$Ch1.Peak.center[i] < unk.part.Meso){ # If they fall below the unknown threshold, assign them "black" color for unknown
fp.data$plastid.color[i] <- 'black'
} else{ #If not unknown, then assign them "blue" color for bluegreen plastid type
if(fp.data$Ch2.Peak.center[i] < bg.part.Meso){
fp.data$plastid.color[i] <- 'blue'
}
}}
dataset.fluo$peak.Ch1 <- fp.data$Ch1.Peak.center[match(dataset.fluo$peak,fp.data$cluster.id)]
dataset.fluo$peak.Ch2 <- fp.data$Ch2.Peak.center[match(dataset.fluo$peak,fp.data$cluster.id)]
dataset.fluo$plastidcolor <- fp.data$plastid.color[match(dataset.fluo$peak,fp.data$cluster.id)]
# Step 5: Map fluorescent data back onto the dataset for export.
dataset.analyzed <- FCDP.Meso.HP
dataset.analyzed$plastidcolor <- dataset.fluo$plastidcolor[match(dataset.analyzed$Capture.ID,dataset.fluo$Capture.ID)]
if(dim(dataset.analyzed[is.na(dataset.analyzed$plastidcolor),])[1]>0){
dataset.analyzed[is.na(dataset.analyzed$plastidcolor),]$plastidcolor <- 'nonfluorescent'}
dataset.analyzed$plastidcolor <- factor(dataset.analyzed$plastidcolor,levels=c('red','blue','black','nonfluorescent'))
# 1. Baseline plot of all data
par(mfrow=c(2,2))
par(mar=c(2,4,2.5,0.5))
plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,ylab='Channel 2 peak (phycoerythrin)')
# 2. Plot of flowpeaks package partitions
par(mar=c(2,2,2.5,2.5))
plot(fp)
# 3. Plot of data with cluster assignments
par(mar=c(4,4,.5,.5))
plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,col=c('red','blue','cyan','green')[as.factor(dataset.fluo$peak)],xlab='Channel 1 peak (chlorophyll)',ylab='Channel 2 peak (phycoerythrin)')
# 4. Plot of data assigned by plastid type
par(mar=c(4,2,.5,2.5))
plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,col=dataset.fluo$plastidcolor,xlab='Channel 1 peak (chlorophyll)')

FCDP.Meso.HP.hasplastid <- dataset.fluo[dataset.fluo$plastidcolor!='black',]
Extract the fluorescent particles – starved
# Step 1: Partition dataset to just fluorescent particles
dataset.fluo <- FCDP.Meso.00[FCDP.Meso.00$Ch1.Peak>0,]
dataset.fluo <- dataset.fluo[dataset.fluo$Ch2.Peak>0,]
# Step 2: Use flowPeaks package to fit and visualize the number of clusters of points
invisible(capture.output(fp <- suppressMessages(flowPeaks(dataset.fluo[,c("Ch1.Peak",'Ch2.Peak')],tol=.1,h0=.4))))
# The "capture.output" function captures the printed output of the flowPeaks function
# The "invisible" function suppresses that output so it doesn't automatically print.
# The "suppressMessages" function suppresses the messages output of flowPeaks
# Step 3: Extract cluster assignments
dataset.fluo$peak <- fp$peaks.cluster
# Step 4: Assign clusters to colors (or unknown)
fp.data <- as.data.frame(summary(fp))
fp.data$plastid.color <- 'red' # Presume plastids are red in color
for(i in 1:dim(fp.data)[1]){
if(fp.data$Ch1.Peak.center[i] < unk.part.Meso){ # If they fall below the unknown threshold, assign them "black" color for unknown
fp.data$plastid.color[i] <- 'black'
} else{ #If not unknown, then assign them "blue" color for bluegreen plastid type
if(fp.data$Ch2.Peak.center[i] < bg.part.Meso){
fp.data$plastid.color[i] <- 'blue'
}
}}
dataset.fluo$peak.Ch1 <- fp.data$Ch1.Peak.center[match(dataset.fluo$peak,fp.data$cluster.id)]
dataset.fluo$peak.Ch2 <- fp.data$Ch2.Peak.center[match(dataset.fluo$peak,fp.data$cluster.id)]
dataset.fluo$plastidcolor <- fp.data$plastid.color[match(dataset.fluo$peak,fp.data$cluster.id)]
# Step 5: Map fluorescent data back onto the dataset for export.
dataset.analyzed <- FCDP.Meso.00
dataset.analyzed$plastidcolor <- dataset.fluo$plastidcolor[match(dataset.analyzed$Capture.ID,dataset.fluo$Capture.ID)]
if(dim(dataset.analyzed[is.na(dataset.analyzed$plastidcolor),])[1]>0){
dataset.analyzed[is.na(dataset.analyzed$plastidcolor),]$plastidcolor <- 'nonfluorescent'}
dataset.analyzed$plastidcolor <- factor(dataset.analyzed$plastidcolor,levels=c('red','blue','black','nonfluorescent'))
# 1. Baseline plot of all data
par(mfrow=c(2,2))
par(mar=c(2,4,2.5,0.5))
plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,ylab='Channel 2 peak (phycoerythrin)')
# 2. Plot of flowpeaks package partitions
par(mar=c(2,2,2.5,2.5))
plot(fp)
# 3. Plot of data with cluster assignments
par(mar=c(4,4,.5,.5))
plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,col=c('red','blue','cyan','green')[as.factor(dataset.fluo$peak)],xlab='Channel 1 peak (chlorophyll)',ylab='Channel 2 peak (phycoerythrin)')
# 4. Plot of data assigned by plastid type
par(mar=c(4,2,.5,2.5))
plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,col=dataset.fluo$plastidcolor,xlab='Channel 1 peak (chlorophyll)')

FCDP.Meso.00.hasplastid <- dataset.fluo[dataset.fluo$plastidcolor!='black',]
xlimits <- c(5.3,6.5)
ylimits <- c(3.5,6.75)
TAto00cols2 <- colorRampPalette(c('indianred3','gray90'))
TAto00cols2 <- colorRampPalette(c(TAcol,'gray80'))
TAto00colset2 <- addalpha(TAto00cols2(length(expt.day)),.5)
TAtoHPcols2 <- colorRampPalette(c('indianred3','gray90','paleturquoise4'))
TAtoHPcols2 <- colorRampPalette(c(TAcol,'gray90',HPcol))
TAtoHPcolset2 <- addalpha(TAtoHPcols2(length(expt.day)),.5)
p1 <- FCDP.Meso.00.hasplastid %>%
ggplot(aes(x=Ch1.Peak, y = Ch2.Peak, color=as.factor(Day)))+
geom_point() +
theme(legend.position='left',
panel.background = element_rect(color='gray20',size=1,fill='white'),
panel.grid.major=element_line(colour="gray80"),
panel.grid.minor=element_line(colour="gray80")) +
scale_color_manual(values=TAto00colset2,name = "Day") +
xlab('Channel 1 peak (chlorophyll fluorescence)') +
ylab('Channel 2 peak (phycoerythrin fluorescence)') +
scale_x_continuous(limits = xlimits) +
scale_y_continuous(limits = ylimits) +
annotate("text", x=5.5, y=6.65, label= expression(paste('(A) Starved')), cex = 5)
#annotate("text", x=6.37, y=3.6, label= expression(paste('(A) Starved')), cex = 5)
p2 <- ggMarginal(p1,type='histogram',groupColour=TRUE,groupFill=TRUE)
p5 <- ggMarginal(p1,groupColour=TRUE)#,groupFill=TRUE)
p3 <- FCDP.Meso.HP.hasplastid %>%
ggplot(aes(x=Ch1.Peak, y = Ch2.Peak, color=as.factor(Day)))+
geom_point() +
theme(legend.position='left',
panel.background = element_rect(color='gray20',size=1,fill='white'),
panel.grid.major=element_line(colour="gray80"),
panel.grid.minor=element_line(colour="gray80")) +
scale_color_manual(values=TAtoHPcolset2,name = "Day") +
xlab('Channel 1 peak (chlorophyll fluorescence)') +
ylab('Channel 2 peak (phycoerythrin fluorescence)') +
scale_x_continuous(limits = xlimits) +
scale_y_continuous(limits = ylimits) +
annotate("text", x=5.63, y=6.65, label= expression(paste('(B) Fed ',italic('H. pacifica'))), cex = 5)
p6 <- ggMarginal(p3,type='histogram',groupColour=TRUE,groupFill=TRUE)
p7 <- ggMarginal(p3,groupColour=TRUE)#,groupFill=TRUE)
plot_grid(p5,p7)

Cell size
FCDP.Meso.00.hasplastid$Volscale <- FCDP.Meso.00.hasplastid$Volume..ABD./1000
FCDP.Meso.HP.hasplastid$Volscale <- FCDP.Meso.HP.hasplastid$Volume..ABD./1000
cellVol.00 <- summarySE(data=FCDP.Meso.00.hasplastid,'Volscale','Day',na.rm=T)
cellVol.HP <- summarySE(data=FCDP.Meso.HP.hasplastid,'Volscale','Day',na.rm=T)
par(mar=c(4,4,1,1))
plot(cellVol.00$Day,cellVol.00$Volscale,type='n',xlab='',ylab='',ylim=c(7,10.500),las=1)
mtext('Experimental Day',side=1,line=2.3)
mtext('Estimated cell volume',side=2,line=2.8)
arrows(cellVol.HP$Day,cellVol.HP$Volscale+cellVol.HP$se,cellVol.HP$Day,cellVol.HP$Volscale-cellVol.HP$se,code=3,angle=90,length=0.01,col='black')
arrows(cellVol.00$Day,cellVol.00$Volscale+cellVol.00$se,cellVol.00$Day,cellVol.00$Volscale-cellVol.00$se,code=3,angle=90,length=0.01,col='black')
lines(cellVol.00$Day,cellVol.00$Volscale)
lines(cellVol.HP$Day,cellVol.HP$Volscale)
points(cellVol.HP$Day,cellVol.HP$Volscale,pch=21,col='black',bg='black',cex=1.5)
points(cellVol.00$Day,cellVol.00$Volscale,pch=21,col='black',bg='white',cex=1.5)
legend(x='topleft',inset=c(0,0),legend=c('Starved',expression(paste('Fed ',italic('H. pacifica')))),bty='n',pch=21,pt.bg=c('white','black'),pt.cex=1.5)

With additional cell properties
par(mar=c(0,4,4,1),mfrow=c(3,1))
FCDP.Meso.00.hasplastid$Volscale <- FCDP.Meso.00.hasplastid$Volume..ABD./1000
FCDP.Meso.HP.hasplastid$Volscale <- FCDP.Meso.HP.hasplastid$Volume..ABD./1000
cellVol.00 <- summarySE(data=FCDP.Meso.00.hasplastid,'Volscale','Day',na.rm=T)
cellVol.HP <- summarySE(data=FCDP.Meso.HP.hasplastid,'Volscale','Day',na.rm=T)
plot(cellVol.00$Day,cellVol.00$Volscale,type='n',xlab='',ylab='',ylim=c(6.5,11),las=1,xaxt='n'); axis(side=1,labels=F)
#mtext('Experimental Day',side=1,line=2.3)
mtext('Estimated cell volume',side=2,line=2.8,cex=.7)
arrows(cellVol.HP$Day,cellVol.HP$Volscale+cellVol.HP$se,cellVol.HP$Day,cellVol.HP$Volscale-cellVol.HP$se,code=3,angle=90,length=0.01,col='black')
arrows(cellVol.00$Day,cellVol.00$Volscale+cellVol.00$se,cellVol.00$Day,cellVol.00$Volscale-cellVol.00$se,code=3,angle=90,length=0.01,col='black')
lines(cellVol.00$Day,cellVol.00$Volscale)
lines(cellVol.HP$Day,cellVol.HP$Volscale)
points(cellVol.HP$Day,cellVol.HP$Volscale,pch=21,col='black',bg='black',cex=1.5)
points(cellVol.00$Day,cellVol.00$Volscale,pch=21,col='black',bg='white',cex=1.5)
legend(x='topleft',legend=c(''),title='(a)',bty='n',cex=1.1)
# Cell perimeter
FCDP.Meso.00.hasplastid$Volscale <- FCDP.Meso.00.hasplastid$Perimeter
FCDP.Meso.HP.hasplastid$Volscale <- FCDP.Meso.HP.hasplastid$Perimeter
cellVol.00 <- summarySE(data=FCDP.Meso.00.hasplastid,'Volscale','Day',na.rm=T)
cellVol.HP <- summarySE(data=FCDP.Meso.HP.hasplastid,'Volscale','Day',na.rm=T)
par(mar=c(2,4,2,1))
plot(cellVol.00$Day,cellVol.00$Volscale,type='n',xlab='',ylab='',las=1,xaxt='n',ylim=c(103,122)); axis(side=1,labels=F)
mtext('Cell perimeter',side=2,line=2.8,cex=.7)
#mtext('Experimental Day',side=1,line=2.3)
arrows(cellVol.HP$Day,cellVol.HP$Volscale+cellVol.HP$se,cellVol.HP$Day,cellVol.HP$Volscale-cellVol.HP$se,code=3,angle=90,length=0.01,col='black')
arrows(cellVol.00$Day,cellVol.00$Volscale+cellVol.00$se,cellVol.00$Day,cellVol.00$Volscale-cellVol.00$se,code=3,angle=90,length=0.01,col='black')
lines(cellVol.00$Day,cellVol.00$Volscale)
lines(cellVol.HP$Day,cellVol.HP$Volscale)
points(cellVol.HP$Day,cellVol.HP$Volscale,pch=21,col='black',bg='black',cex=1.5)
points(cellVol.00$Day,cellVol.00$Volscale,pch=21,col='black',bg='white',cex=1.5)
legend(x='topleft',legend=c(''),title='(b)',bty='n',cex=1.1)
# Length:width ratio
FCDP.Meso.00.hasplastid$Volscale <- FCDP.Meso.00.hasplastid$Length/FCDP.Meso.00.hasplastid$Width
FCDP.Meso.HP.hasplastid$Volscale <- FCDP.Meso.HP.hasplastid$Length/FCDP.Meso.HP.hasplastid$Width
cellVol.00 <- summarySE(data=FCDP.Meso.00.hasplastid,'Volscale','Day',na.rm=T)
cellVol.HP <- summarySE(data=FCDP.Meso.HP.hasplastid,'Volscale','Day',na.rm=T)
par(mar=c(4,4,0,1))
plot(cellVol.00$Day,cellVol.00$Volscale,type='n',xlab='',ylab='',las=1,ylim=c(1.21,1.35))
mtext('Experimental Day',side=1,line=2.3,cex=.7)
mtext('Cell length : width ratio',side=2,line=2.8,cex=.7)
arrows(cellVol.HP$Day,cellVol.HP$Volscale+cellVol.HP$se,cellVol.HP$Day,cellVol.HP$Volscale-cellVol.HP$se,code=3,angle=90,length=0.01,col='black')
arrows(cellVol.00$Day,cellVol.00$Volscale+cellVol.00$se,cellVol.00$Day,cellVol.00$Volscale-cellVol.00$se,code=3,angle=90,length=0.01,col='black')
lines(cellVol.00$Day,cellVol.00$Volscale)
lines(cellVol.HP$Day,cellVol.HP$Volscale)
points(cellVol.HP$Day,cellVol.HP$Volscale,pch=21,col='black',bg='black',cex=1.5)
points(cellVol.00$Day,cellVol.00$Volscale,pch=21,col='black',bg='white',cex=1.5)
legend(x='bottomright',inset=c(0,0),legend=c('Starved',expression(paste('Fed ',italic('H. pacifica')))),bty='n',pch=21,pt.bg=c('white','black'),pt.cex=1.5)
legend(x='topleft',legend=c(''),title='(c)',bty='n',cex=1.1)

Line plots of plastid content
Find composite fluorescence values
ARC$M.meanch1peak <- (ARC$M.pmL.red*ARC$M.Red.meanCh1peak+ARC$M.pmL.blue*ARC$M.Blue.meanCh1peak)/ARC$M.pmL
ARC$M.meanch2peak <- (ARC$M.pmL.red*ARC$M.Red.meanCh2peak+ARC$M.pmL.blue*ARC$M.Blue.meanCh2peak)/ARC$M.pmL
for(i in 1:dim(ARC)[1]){
if(is.na(ARC$M.meanch1peak[i])){
if(is.na(ARC$M.Blue.meanCh1peak[i])){
ARC$M.meanch1peak[i] <- ARC$M.Red.meanCh1peak[i]
}
}
if(is.na(ARC$M.meanch2peak[i])){
if(is.na(ARC$M.Blue.meanCh2peak[i])){
ARC$M.meanch2peak[i] <- ARC$M.Red.meanCh2peak[i]
}
}
if(is.na(ARC$M.meanch1peak[i])){
if(is.na(ARC$M.Red.meanCh1peak[i])){
ARC$M.meanch1peak[i] <- ARC$M.Blue.meanCh1peak[i]
}
}
if(is.na(ARC$M.meanch2peak[i])){
if(is.na(ARC$M.Red.meanCh2peak[i])){
ARC$M.meanch2peak[i] <- ARC$M.Blue.meanCh2peak[i]
}
}
}
Using calibration to turn into prop cell volume with PC and PE
plastids
ARC$prop.PE.predict <- ((ARC$M.meanch2peak-sqrt.PE.b)/sqrt.PE.a)^2
ARC$prop.PC.predict <- (ARC$M.meanch1peak-lm.Chl.a-lm.Chl.c*ARC$prop.PE.predict)/lm.Chl.b
for(i in 1:dim(ARC)[1]){
ARC$prop.PC.predict[i] <- max(0,ARC$prop.PC.predict[i])
}
ARC$prop.plastid.predict <- ARC$prop.PE.predict+ARC$prop.PC.predict
par(mar=c(4,4,1,1),mfrow=c(1,3))
pointcolor <- 'black'
pointcolor.bg <- 'gray80'
plot(jitter(ARC$Day.Numeric,.4),ARC$prop.PE.predict,las=1,xlab='Time (days)',ylab=expression(paste('Prop. cell volume ',italic('T. amphioxeia'),' plastids')),pch=21,bg=c('white',pointcolor.bg)[as.factor(ARC$Prey)],col=pointcolor.bg,cex=1.5)
Summ.Prey <- summarySE(data=ARC, 'prop.PE.predict', groupvars=c('Day.Numeric','Prey'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$prop.PE.predict+Summ.Prey$se,Summ.Prey$Day.Numeric,Summ.Prey$prop.PE.predict-Summ.Prey$se,code=3,angle=90,length=.02)
for(i in 1:length(unique(Summ.Prey$Prey))){
subdat <- Summ.Prey[Summ.Prey$Prey==unique(Summ.Prey$Prey)[i],]
lines(subdat$Day.Numeric,subdat$prop.PE.predict)
}
points(Summ.Prey$Day.Numeric,Summ.Prey$prop.PE.predict,pch=21,cex=1.5,bg=c('white',pointcolor)[as.factor(Summ.Prey$Prey)])
legend(x='bottomleft',inset=c(0.03,0.03),legend=c('Starved',expression(paste('Fed ',italic('H. pacifica')))),pch=21,col='black',pt.bg=c('white',pointcolor),pt.cex=1.5,cex=.8)
text(x=0.5,y=0.6,'(A)')
# Panel B: Phycocyanin plastids
plot(jitter(ARC$Day.Numeric,.4),ARC$prop.PC.predict,las=1,xlab='Time (days)',ylab=expression(paste('Prop. cell volume ',italic('H. pacifica'),' plastids')),pch=21,bg=c('white',pointcolor.bg)[as.factor(ARC$Prey)],col=pointcolor.bg,cex=1.5)
Summ.Prey <- summarySE(data=ARC, 'prop.PC.predict', groupvars=c('Day.Numeric','Prey'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$prop.PC.predict+Summ.Prey$se,Summ.Prey$Day.Numeric,Summ.Prey$prop.PC.predict-Summ.Prey$se,code=3,angle=90,length=.02)
for(i in 1:length(unique(Summ.Prey$Prey))){
subdat <- Summ.Prey[Summ.Prey$Prey==unique(Summ.Prey$Prey)[i],]
lines(subdat$Day.Numeric,subdat$prop.PC.predict)
}
points(Summ.Prey$Day.Numeric,Summ.Prey$prop.PC.predict,pch=21,cex=1.5,bg=c('white',pointcolor)[as.factor(Summ.Prey$Prey)])
text(x=0.5,y=0.207,'(B)')
# Panel C: All Plastids
plot(jitter(ARC$Day.Numeric,.4),ARC$prop.plastid.predict,las=1,xlab='Time (days)',ylab='Prop. cell volume containing plastids',pch=21,bg=c('white',pointcolor.bg)[as.factor(ARC$Prey)],col=pointcolor.bg,cex=1.5,ylim=c(0,0.61))
Summ.Prey <- summarySE(data=ARC, 'prop.plastid.predict', groupvars=c('Day.Numeric','Prey'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$prop.plastid.predict+Summ.Prey$se,Summ.Prey$Day.Numeric,Summ.Prey$prop.plastid.predict-Summ.Prey$se,code=3,angle=90,length=.02)
for(i in 1:length(unique(Summ.Prey$Prey))){
subdat <- Summ.Prey[Summ.Prey$Prey==unique(Summ.Prey$Prey)[i],]
lines(subdat$Day.Numeric,subdat$prop.plastid.predict)
}
points(Summ.Prey$Day.Numeric,Summ.Prey$prop.plastid.predict,pch=21,cex=1.5,bg=c('white',pointcolor)[as.factor(Summ.Prey$Prey)])
text(x=0.5,y=0.6,'(C)')

Population dynamics
# Read in cell count data
dat <- read.csv("/Users/hollyvm/GoogleSync/Mesodinium/Mrubrum_fedHP/LongTermExpt/Mrubrum_fedHP_FunctionalResponseCurves - LTExp_16Dec2023.csv")
# Read in photophysiology data
dat <- dat[!is.na(dat$Flask),]
fire <- read.csv("/Users/hollyvm/GoogleSync/Mesodinium/Mrubrum_fedHP/LongTermExpt/Mrubrum_fedHP_FunctionalResponseCurves - LTExp_PECurves.csv")
fire <- fire[fire$Exp=='LT',]
MRcolumns <- which(substr(colnames(dat),1,6)=='corMR.')
scaleMRcolumns <- which(substr(colnames(dat),1,9)=='altcorMR.')
HPcolumns <- which(substr(colnames(dat),1,3)=='HP.')
Ratiocolumns <- which(substr(colnames(dat),1,6)=='Ratio.')
Chlcolumns <- which(substr(colnames(dat),1,4)=='Chl.')
TotChlcolumns <- which(substr(colnames(dat),1,7)=='TotChl.')
days <- as.numeric(as.character(str_extract(colnames(dat[,MRcolumns]), "\\d+$")))
days[10] <- 8.2
days.chl <- as.numeric(as.character(str_extract(colnames(dat[,Chlcolumns]), "\\d+$")))
treats <- unique(dat$Treat)
flasks <- unique(dat$Flask)
feedpts <- c(2.3,3.3,5.3,7.1,11.1,13.1,14.3,16.3,18.4)
feedvols <- c(0.4,1,2,3,1,3,2,1,1)
dilupts <- c(8)
# Special parameters, etc
Cfix.mult <- 0.01728178 # The multiplier 0.01728178 converts from electrons per chlorophyll molecule per second into pg C per pg chlorophyll-a per hour
par(mar=c(4,4,1,.5),mfrow=c(1,3))
ax.lab.cex <- 0.7
# M. rubrum data
plot(jitter(days,.4),dat[1,MRcolumns],las=1,xlab='',ylab='', type='n',ylim=c(min(dat[,MRcolumns],na.rm=T),max(dat[,MRcolumns],na.rm = T))/1000, xlim=c(min(days)-.1,max(days)+.1),log='y')
abline(v=dilupts,lty=3)
mtext(expression(paste(italic('M. rubrum'),' (',10^3,' cells ',mL^-1,')')),side=2,line=2.2,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(i in 1:length(treats)){
subdat <- dat[dat$Treat==treats[i],]
for(j in 1:dim(dat)[1]){
points(jitter(days,0.4),subdat[j,MRcolumns]/1000,pch=21,bg=TAcols.bg[i],col=TAcol.bg,cex=1.5)
}
lines(days,colMeans(subdat[,MRcolumns],na.rm=T)/1000)
}
for(i in 1:length(days)){
for(j in 1:length(treats)){
subdat <- dat[dat$Treat==treats[j],]
arrows(days[i],mean(subdat[,MRcolumns[i]],na.rm=T)/1000+se(subdat[,MRcolumns[i]])/1000,days[i],mean(subdat[,MRcolumns[i]],na.rm=T)/1000-se(subdat[,MRcolumns[i]])/1000, angle=90, code = 3, length=0.03)
points(days[i],mean(subdat[,MRcolumns[i]],na.rm=T)/1000,pch=21,cex=1.5,col='black',bg=TAcols[j])
}
}
text(dilupts-1.5,.95*rep(min(dat[,MRcolumns]/1000,na.rm=T),length(dilupts)),rep('Dilution',length(dilupts)),pos=4,srt=90)
# H. pacifica data
subdat <- dat[dat$Treat==2,]
plot(jitter(days,.4),subdat[1,HPcolumns],las=1,xlab='',ylab='', type='n',ylim=c(min(subdat[,HPcolumns],na.rm=T),max(subdat[,HPcolumns],na.rm = T))/1000,xlim=c(min(days)-.1,max(days)+.1),log='y')
abline(v=feedpts,lty=2)
abline(v=dilupts,lty=3)
mtext(expression(paste('Cryptophytes (',10^3,' cells ',mL^-1,')')),side=2,line=2,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(j in 1:dim(subdat)[1]){
points(jitter(days,0.4),subdat[j,HPcolumns]/1000,pch=21,bg=HPcols.bg[2],col=HPcol.bg,cex=1.5)
}
lines(days,colMeans(subdat[,HPcolumns],na.rm=T)/1000)
for(i in 1:length(days)){
arrows(days[i],mean(subdat[,HPcolumns[i]],na.rm=T)/1000+se(subdat[,HPcolumns[i]])/1000,days[i],mean(subdat[,HPcolumns[i]],na.rm=T)/1000-se(subdat[,HPcolumns[i]])/1000, angle=90, code = 3, length=0.03)
points(days[i],mean(subdat[,HPcolumns[i]],na.rm=T)/1000,pch=21,cex=1.5,col='black',bg=HPcols[2])
}
text(feedpts+.3,rep(max(subdat[,HPcolumns]/1000)*1.03,length(feedpts)),rep('Pulse fed',length(feedpts)),pos=2,srt=90)
text(dilupts+.0,.93*rep(min(subdat[,HPcolumns]/1000,na.rm=T),length(dilupts)),rep('Dilution',length(dilupts)),pos=4,srt=90)
# Ratio
subdat <- dat[dat$Treat==2,]
plot(jitter(days,.4),subdat[1,Ratiocolumns],las=1,xlab='',ylab='', type='n',ylim=c(min(subdat[,Ratiocolumns],na.rm=T), max(subdat[,Ratiocolumns],na.rm = T)),xlim=c(min(days)-.1,max(days)+.1))
abline(v=feedpts,lty=2)
abline(v=dilupts,lty=3)
mtext(expression(paste('Cryptophytes : ',italic('M. rubrum'))),side=2,line=2.3,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(j in 1:dim(subdat)[1]){
points(jitter(days,0.4),subdat[j,Ratiocolumns],pch=21,bg='gray80',col='gray60',cex=1.5)
}
lines(days,colMeans(subdat[,Ratiocolumns],na.rm=T))
for(i in 1:length(days)){
arrows(days[i],mean(subdat[,Ratiocolumns[i]],na.rm=T)+se(subdat[,Ratiocolumns[i]]),days[i],mean(subdat[,Ratiocolumns[i]],na.rm=T)-se(subdat[,Ratiocolumns[i]]), angle=90, code = 3, length=0.03)
points(days[i],mean(subdat[,Ratiocolumns[i]],na.rm=T),pch=21,cex=1.5,col='black',bg='gray30')
}
text(feedpts+.3,rep(max(subdat[,Ratiocolumns]*1.02,na.rm = T),length(feedpts)),rep('Pulse fed',length(feedpts)),pos=2,srt=90)
text(dilupts+.0,-.3,rep('Dilution',length(dilupts)),pos=4,srt=90)

M. rubrum growth rate
grwindow <- 4 # number of datapoints required for a growth rate calculation
numcalcs <- length(days) - grwindow + 1
gr.data <- as.data.frame(rep(flasks,each=numcalcs))
colnames(gr.data) <- 'Flask'
gr.data$Treat <- dat$Treat[match(gr.data$Flask,dat$Flask)]
gr.data$MeanDay <- NaN
gr.data$GrowthRate <- NaN
gr.data$GrowthRateErr <- NaN
for(i in 1:6){
for(j in 1:numcalcs){
if(sum(days[j:(j+grwindow-1)]%in%c(8,8.2))!=2){ # if the dilution point is not in the string of numbers
gr.data$MeanDay[(i-1)*numcalcs + j] <- mean(days[j:(j+grwindow-1)])
lm1 <- lm(t(log(dat[dat$Flask==flasks[i],MRcolumns[j:(j+grwindow-1)]]))~days[j:(j+grwindow-1)])
gr.data$GrowthRate[(i-1)*numcalcs + j] <- summary(lm1)$coefficients[2,1]
gr.data$GrowthRateErr[(i-1)*numcalcs + j] <- summary(lm1)$coefficients[2,2]
}}}
gr.data <- gr.data[!is.na(gr.data$GrowthRate),]
monochromecols <- c('white','black')
monochromecols.bg <- c('white','gray80')
par(mar=c(4,4,1,1),mfrow=c(1,2))
# M. rubrum data
plot(jitter(days,.4),dat[1,scaleMRcolumns],las=1,xlab='Experimental day',ylab='', type='n',ylim=c(min(dat[,scaleMRcolumns],na.rm=T),max(dat[,scaleMRcolumns],na.rm = T))/1000, xlim=c(min(days)-.1,max(days)+.1),log='y')
#abline(v=dilupts,lty=3)
mtext(expression(paste('Scaled ',italic('M. rubrum'),' popn. (',10^3,' cells ',mL^-1,')')),side=2,line=2,cex=1)
#mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(i in 1:length(treats)){
subdat <- dat[dat$Treat==treats[i],]
for(j in 1:dim(dat)[1]){
points(jitter(days,0.4),subdat[j,scaleMRcolumns]/1000,pch=21,bg=monochromecols.bg[i],col=monochromecols.bg,cex=1.5)
}
lines(days,colMeans(subdat[,scaleMRcolumns],na.rm=T)/1000)
}
for(i in 1:length(days)){
for(j in 1:length(treats)){
subdat <- dat[dat$Treat==treats[j],]
arrows(days[i],mean(subdat[,scaleMRcolumns[i]],na.rm=T)/1000+se(subdat[,scaleMRcolumns[i]])/1000,days[i],mean(subdat[,scaleMRcolumns[i]],na.rm=T)/1000-se(subdat[,scaleMRcolumns[i]])/1000, angle=90, code = 3, length=0.03)
points(days[i],mean(subdat[,scaleMRcolumns[i]],na.rm=T)/1000,pch=21,cex=1.5,col='black',bg=monochromecols[j])
}
}
text(0.25,37,'(A)')
# Growth rates
grwindow <- 7 # number of datapoints required for a growth rate calculation
numcalcs <- length(days) - grwindow + 1
gr.data <- as.data.frame(rep(flasks[1:6],each=numcalcs))
colnames(gr.data) <- 'Flask'
gr.data$Treat <- dat$Treat[match(gr.data$Flask,dat$Flask)]
gr.data$MeanDay <- NaN
gr.data$GrowthRate <- NaN
gr.data$GrowthRateErr <- NaN
for(i in 1:6){
for(j in 1:numcalcs){
gr.data$MeanDay[(i-1)*numcalcs + j] <- mean(days[j:(j+grwindow-1)])
lm1 <- lm(t(log(dat[dat$Flask==flasks[i],scaleMRcolumns[j:(j+grwindow-1)]]))~days[j:(j+grwindow-1)])
gr.data$GrowthRate[(i-1)*numcalcs + j] <- summary(lm1)$coefficients[2,1]
gr.data$GrowthRateErr[(i-1)*numcalcs + j] <- summary(lm1)$coefficients[2,2]
}}
plot(jitter(gr.data$MeanDay,.5),gr.data$GrowthRate,las=1,xlab='Experimental day',ylab=expression(paste(italic('M. rubrum'),' growth rate (',d^-1,')')),pch=21,cex=1.5,bg=monochromecols.bg[as.factor(gr.data$Treat)],col='gray80')
growth.summ <- summarySE(data=gr.data,measurevar='GrowthRate',groupvars=c('MeanDay','Treat'))
#abline(v=dilupts,lty=3)
for(i in 1:length(treats)){
subdat <- growth.summ[growth.summ$Treat==treats[i],]
lines(subdat$MeanDay,subdat$GrowthRate)
arrows(subdat$MeanDay,subdat$GrowthRate + subdat$se, subdat$MeanDay, subdat$GrowthRate - subdat$se, code=3, angle = 90, length = 0.03)
points(subdat$MeanDay,subdat$GrowthRate,pch=21, cex= 1.5, col='black', bg = monochromecols[i])
}
legend(x='bottomleft',inset=c(-0.0,-0.02),legend=c('Starved',expression(paste('Fed ',italic('H. pacifica')))),pch=21,pt.cex=1.5,col='black',pt.bg=monochromecols,bty='n')
text(3.25,0.28,'(B)')

Stacked vertically
gr.data <- gr.data[!is.na(gr.data$GrowthRate),]
monochromecols <- c('white','black')
monochromecols.bg <- c('white','gray80')
par(mar=c(1,4,3,1),mfrow=c(2,1))
# M. rubrum data
plot(jitter(days,.4),dat[1,scaleMRcolumns],las=1,xlab='Experimental day',ylab='', type='n',ylim=c(min(dat[,scaleMRcolumns],na.rm=T),max(dat[,scaleMRcolumns],na.rm = T))/1000, xlim=c(min(days)-.1,max(days)+.1),log='y',xaxt='n')
axis(side=1,labels=F)
#abline(v=dilupts,lty=3)
mtext(expression(paste('Scaled ',italic('M. rubrum'),' popn. (',10^3,' cells ',mL^-1,')')),side=2,line=2,cex=1)
#mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(i in 1:length(treats)){
subdat <- dat[dat$Treat==treats[i],]
for(j in 1:dim(dat)[1]){
points(jitter(days,0.4),subdat[j,scaleMRcolumns]/1000,pch=21,bg=monochromecols.bg[i],col=monochromecols.bg,cex=1.5)
}
lines(days,colMeans(subdat[,scaleMRcolumns],na.rm=T)/1000)
}
for(i in 1:length(days)){
for(j in 1:length(treats)){
subdat <- dat[dat$Treat==treats[j],]
arrows(days[i],mean(subdat[,scaleMRcolumns[i]],na.rm=T)/1000+se(subdat[,scaleMRcolumns[i]])/1000,days[i],mean(subdat[,scaleMRcolumns[i]],na.rm=T)/1000-se(subdat[,scaleMRcolumns[i]])/1000, angle=90, code = 3, length=0.03)
points(days[i],mean(subdat[,scaleMRcolumns[i]],na.rm=T)/1000,pch=21,cex=1.5,col='black',bg=monochromecols[j])
}
}
text(0.25,37,'(A)')
par(mar=c(4,4,0,1))
# Growth rates
grwindow <- 7 # number of datapoints required for a growth rate calculation
numcalcs <- length(days) - grwindow + 1
gr.data <- as.data.frame(rep(flasks[1:6],each=numcalcs))
colnames(gr.data) <- 'Flask'
gr.data$Treat <- dat$Treat[match(gr.data$Flask,dat$Flask)]
gr.data$MeanDay <- NaN
gr.data$GrowthRate <- NaN
gr.data$GrowthRateErr <- NaN
for(i in 1:6){
for(j in 1:numcalcs){
gr.data$MeanDay[(i-1)*numcalcs + j] <- mean(days[j:(j+grwindow-1)])
lm1 <- lm(t(log(dat[dat$Flask==flasks[i],scaleMRcolumns[j:(j+grwindow-1)]]))~days[j:(j+grwindow-1)])
gr.data$GrowthRate[(i-1)*numcalcs + j] <- summary(lm1)$coefficients[2,1]
gr.data$GrowthRateErr[(i-1)*numcalcs + j] <- summary(lm1)$coefficients[2,2]
}}
plot(jitter(gr.data$MeanDay,.5),gr.data$GrowthRate,las=1,xlab='',ylab=expression(paste(italic('M. rubrum'),' growth rate (',d^-1,')')),pch=21,cex=1.5,bg=monochromecols.bg[as.factor(gr.data$Treat)],col='gray80')
mtext('Experimental day',side=1,line=2.3)
growth.summ <- summarySE(data=gr.data,measurevar='GrowthRate',groupvars=c('MeanDay','Treat'))
#abline(v=dilupts,lty=3)
for(i in 1:length(treats)){
subdat <- growth.summ[growth.summ$Treat==treats[i],]
lines(subdat$MeanDay,subdat$GrowthRate)
arrows(subdat$MeanDay,subdat$GrowthRate + subdat$se, subdat$MeanDay, subdat$GrowthRate - subdat$se, code=3, angle = 90, length = 0.03)
points(subdat$MeanDay,subdat$GrowthRate,pch=21, cex= 1.5, col='black', bg = monochromecols[i])
}
legend(x='bottomleft',inset=c(-0.0,-0.02),legend=c('Starved',expression(paste('Fed ',italic('H. pacifica')))),pch=21,pt.cex=1.5,col='black',pt.bg=monochromecols,bty='n')
text(3.25,0.28,'(B)')

Photophysiology
Fit PE curves to all data
flasks <- unique(fire$Flask)
days.fire <- unique(fire$ExpDay)
fire$Flask.Day <- paste(fire$Flask,fire$ExpDay,sep='.')
# Create PE dataframe
PEres <- as.data.frame(cbind(rep(flasks,length(days.fire)),rep(days.fire,each=length(flasks))))
colnames(PEres) <- c('Flask','ExpDay')
PEres$Flask.Day <- paste(PEres$Flask,PEres$ExpDay,sep='.')
PEres$Treatment <- fire$Treatment[match(PEres$Flask,fire$Flask)]
PEres$Prey <- fire$Prey[match(PEres$Flask,fire$Flask)]
PEres$Rep <- fire$Rep[match(PEres$Flask,fire$Flask)]
PEres$TargetHP <- fire$TargetHP[match(PEres$Flask,fire$Flask)]
PEres$FvFm <- fire[fire$PAR==0,]$Fv.Fm[match(PEres$Flask.Day,fire[fire$PAR==0,]$Flask.Day)]
PEres$P.est <- NA
PEres$a.est <- NA
PEres$Pmax <- NA
PEres$Pmax.err <- NA
PEres$Pmax.p <- NA
PEres$a <- NA
PEres$a.err <- NA
PEres$a.p <- NA
PEres$Ek <- NA
PEres$P_I <- NA
PEres$P_I.C <- NA
PEres$P.interp <- NA
PARset <- seq(from =0, to = max(fire$PAR,na.rm=T), length.out = 100) # Create a holding vector of light levels to use for plotting later
# Fit curves
for(i in 1:(dim(PEres)[1])){
if((i-1)%%30 == 0){ # plot results in 30-panel figures
par(mar=c(1,1,1,1),mfrow=c(5,6))
}
subdat <- fire[fire$Flask.Day==PEres$Flask.Day[i],] # Subset the data
if(!is.na(subdat$PAR[1])){
plot(subdat$PAR, subdat$ETR,xlim=c(0,300),ylim=c(0,300)); text(max(PARset)*0.5,5,PEres$Flask.Day[i]) # Plot the raw data
PEres$FvFm[i] <- subdat[subdat$PAR==0,]$Fv.Fm # Save FvFm
for(j in 1:dim(subdat)[1]){if(subdat$ETR[j]<0){subdat$ETR[j]<- NaN}} # Remove negative values for ETR
# note: usually points with negative qpcoef are those that are difficult to fit. in some cases, you can comment the following line of code to include those points
for(j in 1:dim(subdat)[1]){if(subdat$qp.coeff[j]<0){subdat$ETR[j]<- NaN}} # Remove negative values for qp coef
#for(j in 5:dim(subdat[1])){if(!is.nan(subdat$ETR[j-1])){if(!is.nan(subdat$ETR[j])){
# if(subdat$ETR[j] > 1.3*subdat$ETR[j-1]){subdat$ETR[j]<-NaN}}}
P.est <- max(subdat$ETR,na.rm=TRUE)
PEres$P.est[i] <- P.est
a.est <- lm(subdat$ETR[1:3]~subdat$PAR[1:3])$coef[2]
PEres$a.est[i] <- a.est
points(subdat$PAR,subdat$ETR,col='red')
rm(nls.summ)
try({
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= subdat, start = list(P = P.est, a = a.est),nls.control(tol=1e-7))) }) # Do curve fit; 'try' wrapper prevents an error from causing the program to crash
if(exists("nls.summ")){ # Only save these data if the curve fit was successful
PEres$Pmax[i] <- nls.summ$p[1,1] # Save various statistical parameters
PEres$Pmax.err[i] <- nls.summ$p[1,2]
PEres$Pmax.p[i] <- nls.summ$p[1,4]
PEres$a[i] <- nls.summ$p[2,1]
PEres$a.err[i] <- nls.summ$p[2,2]
PEres$a.p[i] <- nls.summ$p[2,4]
PEres$P.interp[i] <- interp1(subdat$PAR,subdat$ETR,xi=20)
predictset <- PEres$Pmax[i]*tanh(PEres$a[i]*PARset/PEres$Pmax[i])
lines(PARset,predictset, col='royalblue3') # Overlay the curve fit
predictset2 <- P.est*tanh(a.est*PARset/P.est)
lines(PARset,predictset2, col='green') # Overlay the curve fit
}}}
## Warning in rm(nls.summ): object 'nls.summ' not found


PEres$Ek <- PEres$Pmax / PEres$a
PEres$P_I <- PEres$Pmax * tanh(PEres$a * 20 / PEres$Pmax)
PEres$P_I.C <- PEres$P_I* Cfix.mult
PEres$P.interp.C <- PEres$P.interp* Cfix.mult
# Export data
# # # write.csv(PEres,'/Users/hollyvm/GoogleSync/Mesodinium/Mrubrum_fedHP/LongTermExpt/PEres.csv')
We’ll incorporate data from reference cryptophyte PE curves as well
HP.ref <- read.csv("/Users/hollyvm/GoogleSync/Mesodinium/Mrubrum_fedHP/Mrubrum_fedHP_FunctionalResponseCurves - Exp2_PECurves.csv")
HP.ref <- HP.ref[HP.ref$TargetMR==0&HP.ref$Treatment%in%c(3,5),]
HP.ref <- read.csv("/Users/hollyvm/GoogleSync/Mesodinium/Mrubrum_fedHP/LongTermExpt/Mrubrum_fedHP_FunctionalResponseCurves - LTExp_PECurves.csv")
TA.ref <- HP.ref[HP.ref$Treatment=='TA',]
HP.ref <- HP.ref[HP.ref$Treatment=='HP',]
# Color spectrum for experimental days
TAto00cols2 <- colorRampPalette(c(TAcol,'gray80'))
TAto00colset2 <- addalpha(TAto00cols2(length(days.fire)),.5)
TAto00colset <- TAto00cols2(length(days.fire))
TAtoHPcols2 <- colorRampPalette(c(TAcol,'gray90',HPcol))
TAtoHPcolset2 <- addalpha(TAtoHPcols2(length(expt.day)),.5)
TAtoHPcolset <- TAtoHPcols2(length(expt.day))
PARset2 <- seq(from = 0, to = 340, length.out=200)
par(mar=c(4,3.5,1,0),mfrow=c(1,2))
# Panel 1: Treatment 1
subdat <- PEres[PEres$Treatment==1,]
subdat.fire <- fire[fire$Treatment==1,]
plot(subdat.fire$PAR,subdat.fire$ETR,xlim=c(0,333),ylim=c(0,295),type='n',xlab='',ylab='',las=1)
mtext(expression(paste('Electron transport rate (e- chl-',italic(a)^-1,' ',s^-1,')')),side=2,line=2.5)
for(i in 1:length(days.fire)){
subdat.fire <- fire[fire$Treatment==1 & fire$ExpDay==days.fire[i],]
points(subdat.fire$PAR,subdat.fire$ETR,pch=21,bg=TAto00colset2[i],col=TAto00colset2[i])
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= subdat.fire, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset/Pmax.fit)
lines(PARset,predictset, col=TAto00colset[i],lwd=2) # Overlay the curve fit
}
text(x=-10,y=285,'(A) Starved',pos=4)
text(x=293,y=265,'Day 0',pos=4,cex=.7)
text(x=293,y=193,'Day 10',pos=4,cex=.7,col=TAto00colset[7])
text(x=293,y=93,'Day 22',pos=4,cex=.7,col=TAto00colset[length(days.fire)])
# Add reference line from HP data
#points(HP.ref$PAR,HP.ref$ETR,col=addalpha(HPcol.bg,.2),bg=addalpha(HPcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= HP.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, col='gray30',lwd=2,lty=3)
text(335,Pmax.fit*.88,expression(paste('Free-living ',italic('H. pacifica'))),pos=2,cex=.7,col='black')
# Add reference line from TA data
#points(TA.ref$PAR,TA.ref$ETR,col=addalpha(TAcol.bg,.2),bg=addalpha(TAcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= TA.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, col='gray30',lwd=2,lty=2)
text(335,Pmax.fit*.96,expression(paste('Free-living ',italic('T. amphioxeia'))),pos=2,cex=.7,col='black')
# Panel 1: Treatment 1
par(mar=c(4,1,1,2.5))
subdat.fire <- fire[fire$Treatment==2,]
plot(subdat.fire$PAR,subdat.fire$ETR,xlim=c(0,333),ylim=c(0,295),type='n',xlab='',ylab='',las=1,yaxt='n')
axis(side=2,labels=F)
mtext(expression(paste('Light level (μmol quanta ',m^-2,' ',s^-1,')')),side=1,line=2.3,at=-30)
for(i in 1:length(days.fire)){
subdat.fire <- fire[fire$Treatment==2 & fire$ExpDay==days.fire[i],]
points(subdat.fire$PAR,subdat.fire$ETR,pch=21,bg=TAtoHPcolset2[i],col=TAtoHPcolset2[i])
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= subdat.fire, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset/Pmax.fit)
lines(PARset,predictset, col=TAtoHPcolset[i],lwd=2) # Overlay the curve fit
}
text(x=-10,y=285,expression(paste('(B) Fed ',italic('H. pacifica'))),pos=4)
text(x=293,y=252,'Day 0',pos=4,cex=.7)
text(x=293,y=170,'Day 10',pos=4,cex=.7,col=TAto00colset[7])
text(x=293,y=60,'Day 22',pos=4,cex=.7,col=TAto00colset[length(days.fire)])
# Add reference line from HP data
#points(HP.ref$PAR,HP.ref$ETR,col=addalpha(HPcol.bg,.2),bg=addalpha(HPcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= HP.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, lwd=2,lty=3, col='gray30') #col=addalpha(HPcol.bg,.8),
#text(315,Pmax.fit*.85,expression(paste('Free-living ',italic('H. pacifica'))),pos=2,cex=.85,col='black')
HPfit <- nls.summ
# Add reference line from TA data
#points(TA.ref$PAR,TA.ref$ETR,col=addalpha(TAcol.bg,.2),bg=addalpha(TAcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= TA.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, col='gray30',lwd=2,lty=2)

#text(315,Pmax.fit*.93,expression(paste('Free-living ',italic('T. amphioxeia'))),pos=2,cex=.85,col='black')
TAfit <- nls.summ
# Color spectrum for experimental days
TAto00cols2 <- colorRampPalette(c('black','gray80'))
TAto00colset2 <- addalpha(TAto00cols2(length(days.fire)),.5)
TAto00colset <- TAto00cols2(length(days.fire))
TAtoHPcols2 <- colorRampPalette(c('black','gray90'))
TAtoHPcolset2 <- addalpha(TAtoHPcols2(length(expt.day)),.5)
TAtoHPcolset <- TAtoHPcols2(length(expt.day))
PARset2 <- seq(from = 0, to = 340, length.out=200)
par(mar=c(4,3.5,1,0),mfrow=c(1,2))
# Panel 1: Treatment 1
subdat <- PEres[PEres$Treatment==1,]
subdat.fire <- fire[fire$Treatment==1,]
plot(subdat.fire$PAR,subdat.fire$ETR,xlim=c(0,333),ylim=c(0,295),type='n',xlab='',ylab='',las=1)
mtext(expression(paste('Electron transport rate (e- chl-',italic(a)^-1,' ',s^-1,')')),side=2,line=2.5)
for(i in 1:length(days.fire)){
subdat.fire <- fire[fire$Treatment==1 & fire$ExpDay==days.fire[i],]
points(subdat.fire$PAR,subdat.fire$ETR,pch=21,bg=TAto00colset2[i],col=TAto00colset2[i])
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= subdat.fire, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset/Pmax.fit)
lines(PARset,predictset, col=TAto00colset[i],lwd=2) # Overlay the curve fit
}
text(x=-10,y=285,'Starved',pos=4)
text(x=293,y=265,'Day 0',pos=4,cex=.7)
text(x=293,y=193,'Day 10',pos=4,cex=.7,col=TAto00colset[7])
text(x=293,y=93,'Day 22',pos=4,cex=.7,col=TAto00colset[length(days.fire)])
# Add reference line from HP data
#points(HP.ref$PAR,HP.ref$ETR,col=addalpha(HPcol.bg,.2),bg=addalpha(HPcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= HP.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, col='gray30',lwd=2,lty=2)
text(335,Pmax.fit*.88,expression(paste('Free-living ',italic('H. pacifica'))),pos=2,cex=.7,col='black')
# Add reference line from TA data
#points(TA.ref$PAR,TA.ref$ETR,col=addalpha(TAcol.bg,.2),bg=addalpha(TAcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= TA.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, col='gray30',lwd=2,lty=2)
text(335,Pmax.fit*.96,expression(paste('Free-living ',italic('T. amphioxeia'))),pos=2,cex=.7,col='black')
# Panel 1: Treatment 1
par(mar=c(4,1,1,2.5))
subdat.fire <- fire[fire$Treatment==2,]
plot(subdat.fire$PAR,subdat.fire$ETR,xlim=c(0,333),ylim=c(0,295),type='n',xlab='',ylab='',las=1,yaxt='n')
axis(side=2,labels=F)
mtext(expression(paste('Light level (μmol quanta ',m^-2,' ',s^-1,')')),side=1,line=2.3,at=-30)
for(i in 1:length(days.fire)){
subdat.fire <- fire[fire$Treatment==2 & fire$ExpDay==days.fire[i],]
points(subdat.fire$PAR,subdat.fire$ETR,pch=21,bg=TAtoHPcolset2[i],col=TAtoHPcolset2[i])
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= subdat.fire, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset/Pmax.fit)
lines(PARset,predictset, col=TAtoHPcolset[i],lwd=2) # Overlay the curve fit
}
text(x=-10,y=285,expression(paste('Fed ',italic('H. pacifica'))),pos=4)
text(x=293,y=252,'Day 0',pos=4,cex=.7)
text(x=293,y=170,'Day 10',pos=4,cex=.7,col=TAto00colset[7])
text(x=293,y=60,'Day 22',pos=4,cex=.7,col=TAto00colset[length(days.fire)])
# Add reference line from HP data
#points(HP.ref$PAR,HP.ref$ETR,col=addalpha(HPcol.bg,.2),bg=addalpha(HPcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= HP.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, lwd=2,lty=2, col='gray30') #col=addalpha(HPcol.bg,.8),
#text(315,Pmax.fit*.85,expression(paste('Free-living ',italic('H. pacifica'))),pos=2,cex=.85,col='black')
HPfit <- nls.summ
# Add reference line from TA data
#points(TA.ref$PAR,TA.ref$ETR,col=addalpha(TAcol.bg,.2),bg=addalpha(TAcol.bg,.2),pch=21)
nls.summ <- summary(nls(ETR ~ P * tanh(a * PAR / P), data= TA.ref, start = list(P = P.est, a = a.est),nls.control(tol=1e-7)))
Pmax.fit <- nls.summ$p[1,1] # Save various statistical parameters
a.fit <- nls.summ$p[2,1]
predictset <- Pmax.fit*tanh(a.fit*PARset2/Pmax.fit)
lines(PARset2,predictset, col='gray30',lwd=2,lty=2)

#text(315,Pmax.fit*.93,expression(paste('Free-living ',italic('T. amphioxeia'))),pos=2,cex=.85,col='black')
TAfit <- nls.summ
Chlorophyll content and photophysiology metrics.
Chlcols <- c('white','black')
Chlcols.bg <- c('white','gray70')
Chlcol.bg <- 'gray70'
par(mfrow=c(2,3))
ax.lab.cex <- 0.7
jittamt <- 1
# First row
par(mar=c(1,4,4,.5))
# Chlorophyll per cell
plot(jitter(days.chl,jittamt),dat[1,Chlcolumns],las=1,xlab='',ylab='', type='n',ylim=c(min(dat[,Chlcolumns],na.rm=T),max(dat[,Chlcolumns],na.rm = T)), xlim=c(min(days.chl)-.1,max(days.chl)+.1),xaxt='n'); axis(side=1,labels=F)
mtext(expression(paste('Cellular chlorophyll (pg chl-',italic('a'),' ',cell^-1,')')),side=2,line=2.2,cex=ax.lab.cex)
#mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(i in 1:length(treats)){
subdat <- dat[dat$Treat==treats[i],]
for(j in 1:dim(dat)[1]){
points(jitter(days.chl,0.4),subdat[j,Chlcolumns],pch=21,bg=Chlcols.bg[i],col=Chlcol.bg,cex=1.5)
}
lines(days.chl,colMeans(subdat[,Chlcolumns]))
}
for(i in 1:length(days.chl)){
for(j in 1:length(treats)){
subdat <- dat[dat$Treat==treats[j],]
arrows(days.chl[i],mean(subdat[,Chlcolumns[i]])+se(subdat[,Chlcolumns[i]]),days.chl[i],mean(subdat[,Chlcolumns[i]])-se(subdat[,Chlcolumns[i]]), angle=90, code = 3, length=0.03)
points(days.chl[i],mean(subdat[,Chlcolumns[i]]),pch=21,cex=1.5,col='black',bg=Chlcols[j])
}
}
legend(x='topright',inset=c(.01,.01), legend=c('Starved',expression(paste('Fed ',italic('H. pacifica')))),pch=21,pt.bg=Chlcols,pt.cex=1.5, bty='n')
text(0.3,33,'(A)')
# # Chlorophyll per mL
# Chlcol.bg <- 'black'
#
# plot(jitter(days.chl,jittamt),dat[1,TotChlcolumns],las=1,xlab='',ylab='', type='n',ylim=c(min(dat[,TotChlcolumns],na.rm=T),max(dat[,TotChlcolumns],na.rm = T)), xlim=c(min(days.chl)-.1,max(days.chl)+.1))
# mtext(expression(paste('Total chlorophyll (ng chl-',italic('a'),' ',mL^-1,')')),side=2,line=2.2,cex=ax.lab.cex)
# mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
# for(i in 1:length(treats)){
# subdat <- dat[dat$Treat==treats[i],]
# for(j in 1:dim(dat)[1]){
# points(jitter(days.chl,0.4),subdat[j,TotChlcolumns],pch=21,bg=Chlcols.bg[i],col=Chlcol.bg,cex=1.5)
# }
# lines(days.chl,colMeans(subdat[,TotChlcolumns]))
# }
# for(i in 1:length(days.chl)){
# for(j in 1:length(treats)){
# subdat <- dat[dat$Treat==treats[j],]
# arrows(days.chl[i],mean(subdat[,TotChlcolumns[i]])+se(subdat[,TotChlcolumns[i]]),days.chl[i],mean(subdat[,TotChlcolumns[i]])-se(subdat[,TotChlcolumns[i]]), angle=90, code = 3, length=0.03)
# points(days.chl[i],mean(subdat[,TotChlcolumns[i]]),pch=21,cex=1.5,col='black',bg=Chlcols[j])
# }
# }
# Fv/Fm
plot(jitter(PEres$ExpDay,jittamt),PEres$FvFm,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1),ylim=c(.41,.58),xaxt='n'); axis(side=1,labels=F)
mtext(expression(paste('Photosynthetic efficiency, (',F[v],' / ',F[m],')')),side=2,line=2.6,cex=ax.lab.cex)
#mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'FvFm',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$FvFm)
arrows(subdat$ExpDay, subdat$FvFm+subdat$se, subdat$ExpDay, subdat$FvFm-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$FvFm, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# HP reference
abline(h = mean(HP.ref[HP.ref$PAR==0,]$Fv.Fm),lty=3)
# TA reference
abline(h = mean(TA.ref[HP.ref$PAR==0,]$Fv.Fm),lty=2)
points(0.3,.58,pch=22,bg='white',col='white',cex=2); text(0.3,0.575,'(B)')
# P_I.C
plot(jitter(PEres$ExpDay,jittamt),PEres$P_I.C,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1),ylim=c(.3,.58),xaxt='n')
axis(side=1,labels=F)
mtext(expression(paste('Phot. rate, ',P[I][,][C],' (pg C pg chl-',italic('a')^-1,' h',r^-1,')')),side=2,line=2.6,cex=ax.lab.cex)
#mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'P_I.C',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$P_I.C)
arrows(subdat$ExpDay, subdat$P_I.C+subdat$se, subdat$ExpDay, subdat$P_I.C-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$P_I.C, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# HP reference
abline(h = HPfit$p[1,1] * tanh(HPfit$p[2,1] * 20 / HPfit$p[1,1])* Cfix.mult,lty=3)
# TA reference
abline(h = TAfit$p[1,1] * tanh(TAfit$p[2,1] * 20 / TAfit$p[1,1])* Cfix.mult,lty=2)
points(0.3,.57,pch=21,bg='white',col='white',cex=3); text(0.3,0.57,'(C)')
# Second row
par(mar=c(4,4,1,.5))
# Pmax
#Chlcol.bg <- 'gray70'
plot(jitter(PEres$ExpDay,jittamt),PEres$Pmax,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1))
mtext(expression(paste('Max. phot., ',P[m][a][x],' (e- molec. chl-',italic('a')^-1,' ',s^-1,')')),side=2,line=2.6,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'Pmax',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$Pmax)
arrows(subdat$ExpDay, subdat$Pmax+subdat$se, subdat$ExpDay, subdat$Pmax-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$Pmax, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# HP reference
abline(h = HPfit$p[1,1],lty=3)
# TA reference
abline(h = TAfit$p[1,1],lty=2)
points(0.3,293,pch=21,bg='white',col='white',cex=3); text(0.3,293,'(D)')
# a
plot(jitter(PEres$ExpDay,jittamt),PEres$a,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1),ylim=c(0.9,1.8))
mtext(expression(paste('Light sensitivity, a')),side=2,line=2.9,cex=ax.lab.cex)
mtext(expression(paste('(e- ',m^2,' μmol chl-',italic('a')^-1,' ',quanta^-1,')')),side=2,line=1.9,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'a',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$a)
arrows(subdat$ExpDay, subdat$a+subdat$se, subdat$ExpDay, subdat$a-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$a, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# HP reference
abline(h = HPfit$p[2,1],lty=3)
# TA reference
abline(h = TAfit$p[2,1],lty=2)
points(0.3,1.785,pch=21,bg='white',col='white',cex=3); text(0.3,1.785,'(E)')
# E_k
plot(jitter(PEres$ExpDay,jittamt),PEres$Ek,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1))
mtext(expression(paste('Saturating light intensity, ',E[k])),side=2,line=2.6,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'Ek',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$Ek)
arrows(subdat$ExpDay, subdat$Ek+subdat$se, subdat$ExpDay, subdat$Ek-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$Ek, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# HP reference
abline(h = HPfit$p[1,1]/HPfit$p[2,1],lty=3)
# TA reference
abline(h = TAfit$p[1,1]/TAfit$p[2,1],lty=2)
text(-1.5,HPfit$p[1,1]/HPfit$p[2,1]*1.07,expression(paste('Free-living ',italic('H. pacifica'))),pos=4,cex=1,col='black')
text(23.5,TAfit$p[1,1]/TAfit$p[2,1]*1.02,expression(paste('Free-living ',italic('T. amphioxeia'))),pos=2,cex=1,col='black')
points(0.3,192,pch=21,bg='white',col='white',cex=3); text(0.3,192,'(F)')

Cleaner version
Chlcols <- c('white','black')
Chlcols.bg <- c('white','gray70')
Chlcol.bg <- 'gray70'
par(mar=c(4,4,1,.5),mfrow=c(2,3))
ax.lab.cex <- 0.7
jittamt <- 1
# Chlorophyll per cell
plot(jitter(days.chl,jittamt),dat[1,Chlcolumns],las=1,xlab='',ylab='', type='n', xlim=c(min(days.chl)-.1,max(days.chl)+.1), ylim=c(10,32))
mtext(expression(paste('Cellular chlorophyll (pg chl-',italic('a'),' ',cell^-1,')')),side=2,line=2.2,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(i in 1:length(treats)){
subdat <- dat[dat$Treat==treats[i],]
for(j in 1:dim(dat)[1]){
#points(jitter(days.chl,0.4),subdat[j,Chlcolumns],pch=21,bg=Chlcols.bg[i],col=Chlcol.bg,cex=1.5)
}
lines(days.chl,colMeans(subdat[,Chlcolumns]))
}
for(i in 1:length(days.chl)){
for(j in 1:length(treats)){
subdat <- dat[dat$Treat==treats[j],]
arrows(days.chl[i],mean(subdat[,Chlcolumns[i]])+se(subdat[,Chlcolumns[i]]),days.chl[i],mean(subdat[,Chlcolumns[i]])-se(subdat[,Chlcolumns[i]]), angle=90, code = 3, length=0.03)
points(days.chl[i],mean(subdat[,Chlcolumns[i]]),pch=21,cex=1.5,col='black',bg=Chlcols[j])
}
}
legend(x='topright',inset=c(.01,.01), legend=c('Starved',expression(paste('Fed ',italic('H. pacifica')))),pch=21,pt.bg=Chlcols,pt.cex=1.5)
# Chlorophyll per mL
plot(jitter(days.chl,jittamt),dat[1,TotChlcolumns],las=1,xlab='',ylab='', type='n',ylim=c(min(dat[,TotChlcolumns],na.rm=T),max(dat[,TotChlcolumns],na.rm = T)), xlim=c(min(days.chl)-.1,max(days.chl)+.1))
mtext(expression(paste('Total chlorophyll (ng chl-',italic('a'),' ',mL^-1,')')),side=2,line=2.4,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
for(i in 1:length(treats)){
subdat <- dat[dat$Treat==treats[i],]
for(j in 1:dim(dat)[1]){
#points(jitter(days.chl,0.4),subdat[j,TotChlcolumns],pch=21,bg=Chlcols.bg[i],col=Chlcol.bg,cex=1.5)
}
lines(days.chl,colMeans(subdat[,TotChlcolumns]))
}
for(i in 1:length(days.chl)){
for(j in 1:length(treats)){
subdat <- dat[dat$Treat==treats[j],]
arrows(days.chl[i],mean(subdat[,TotChlcolumns[i]])+se(subdat[,TotChlcolumns[i]]),days.chl[i],mean(subdat[,TotChlcolumns[i]])-se(subdat[,TotChlcolumns[i]]), angle=90, code = 3, length=0.03)
points(days.chl[i],mean(subdat[,TotChlcolumns[i]]),pch=21,cex=1.5,col='black',bg=Chlcols[j])
}
}
# Fv/Fm
plot(jitter(PEres$ExpDay,jittamt),PEres$FvFm,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1),type='n',ylim=c(0.42,0.58))
# HP reference
abline(h = mean(HP.ref[HP.ref$PAR==0,]$Fv.Fm),lty=3)
# TA reference
abline(h = mean(TA.ref[HP.ref$PAR==0,]$Fv.Fm),lty=2)
mtext(expression(paste('Photosynthetic efficiency, (',F[v],' / ',F[m],')')),side=2,line=2.6,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'FvFm',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$FvFm)
arrows(subdat$ExpDay, subdat$FvFm+subdat$se, subdat$ExpDay, subdat$FvFm-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$FvFm, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# P_I.C
plot(jitter(PEres$ExpDay,jittamt),PEres$P_I.C,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1),type='n',ylim=c(.32,.58))
# HP reference
abline(h = HPfit$p[1,1] * tanh(HPfit$p[2,1] * 20 / HPfit$p[1,1])* Cfix.mult,lty=3)
# TA reference
abline(h = TAfit$p[1,1] * tanh(TAfit$p[2,1] * 20 / TAfit$p[1,1])* Cfix.mult,lty=2)
mtext(expression(paste('Phot. rate, ',P[I][,][C],' (pg C pg chl-',italic('a')^-1,' h',r^-1,')')),side=2,line=2.6,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'P_I.C',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$P_I.C)
arrows(subdat$ExpDay, subdat$P_I.C+subdat$se, subdat$ExpDay, subdat$P_I.C-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$P_I.C, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# Pmax
plot(jitter(PEres$ExpDay,jittamt),PEres$Pmax,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1),type='n')
# HP reference
abline(h = HPfit$p[1,1],lty=3)
# TA reference
abline(h = TAfit$p[1,1],lty=2)
mtext(expression(paste('Max. phot., ',P[m][a][x],' (e- molec. chl-',italic('a')^-1,' ',s^-1,')')),side=2,line=2.4,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
FvFmsumm <- summarySE(data = PEres, 'Pmax',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$Pmax)
arrows(subdat$ExpDay, subdat$Pmax+subdat$se, subdat$ExpDay, subdat$Pmax-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$Pmax, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
# # a
# plot(jitter(PEres$ExpDay,jittamt),PEres$a,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1))
# mtext(expression(paste('Light sensitivity, a')),side=2,line=2.6,cex=ax.lab.cex)
# mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
# FvFmsumm <- summarySE(data = PEres, 'a',groupvars=c('ExpDay','Treatment'))
# for(i in 1:length(treats)){
# subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
# lines(subdat$ExpDay,subdat$a)
# arrows(subdat$ExpDay, subdat$a+subdat$se, subdat$ExpDay, subdat$a-subdat$se, angle=90, code = 3, length=0.03)
# points(subdat$ExpDay,subdat$a, cex=1.5, col='black',bg=Chlcols[i], pch=21)
# }
# E_k
plot(jitter(PEres$ExpDay,jittamt),PEres$Ek,las=1,xlab='',ylab='', pch=21, col=Chlcol.bg, bg = Chlcols.bg[as.factor(PEres$Treatment)], cex=1.5, xlim=c(min(PEres$ExpDay)-.1,max(PEres$ExpDay)+.1),type='n')
mtext(expression(paste('Saturating light intensity, ',E[k])),side=2,line=2.6,cex=ax.lab.cex)
mtext('Time (days)',side=1,line=2,cex=ax.lab.cex)
# HP reference
abline(h = HPfit$p[1,1]/HPfit$p[2,1],lty=3)
# TA reference
abline(h = TAfit$p[1,1]/TAfit$p[2,1],lty=2)
text(-1.5,HPfit$p[1,1]/HPfit$p[2,1]*1.07,expression(paste('Free-living ',italic('H. pacifica'))),pos=4,cex=1,col='black')
text(23.5,TAfit$p[1,1]/TAfit$p[2,1]*1.02,expression(paste('Free-living ',italic('T. amphioxeia'))),pos=2,cex=1,col='black')
FvFmsumm <- summarySE(data = PEres, 'Ek',groupvars=c('ExpDay','Treatment'))
for(i in 1:length(treats)){
subdat <- FvFmsumm[FvFmsumm$Treatment==treats[i],]
lines(subdat$ExpDay,subdat$Ek)
arrows(subdat$ExpDay, subdat$Ek+subdat$se, subdat$ExpDay, subdat$Ek-subdat$se, angle=90, code = 3, length=0.03)
points(subdat$ExpDay,subdat$Ek, cex=1.5, col='black',bg=Chlcols[i], pch=21)
}
